Chapter 8 HMSC analysis

8.1 Load data

load("data/data.Rdata")
load("hmsc/fit_model1_250_10.Rdata")

8.2 Variance partitioning

# Compute variance partitioning
varpart=computeVariancePartitioning(m)

varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
   group_by(variable) %>%
   summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
   tt()
tinytable_tt9i5ahgkx9jhlpi5kqy
variable mean sd
Random: location 37.900015 25.317903
logseqdepth 56.110626 25.796874
sex 4.937460 5.612719
origin 1.051899 1.282563
# Basal tree
varpart_tree <- genome_tree

#Varpart table
varpart_table <- varpart$vals %>%
   as.data.frame() %>%
   rownames_to_column(var="variable") %>%
   pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
   mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
   mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)

colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__"))%>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% varpart_tree$tip.label) %>%
    arrange(match(genome, varpart_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
varpart_tree <- varpart_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
   scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()

# Add variance stacked barplot
vertical_tree <-  varpart_tree +
       scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
        geom_fruit(
             data=varpart_table,
             geom=geom_bar,
             mapping = aes(x=value, y=genome, fill=variable, group=variable),
             pwidth = 2,
             offset = 0.05,
             width= 1,
             orientation="y",
             stat="identity")+
      labs(fill="Variable")

vertical_tree

8.3 Posterior estimates

# Select desired support threshold
support=0.9
negsupport=1-support

# Basal tree
postestimates_tree <- genome_tree

# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
    mutate(value = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
    pivot_wider(names_from = variable, values_from = value) %>%
    #select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
    column_to_rownames(var="genome")

#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
    mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
    column_to_rownames(var = "genome") %>%
    dplyr::select(phylum)


colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
    mutate(phylum=str_remove_all(phylum, "p__")) %>%
    right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
    filter(genome %in% postestimates_tree$tip.label) %>%
    arrange(match(genome, postestimates_tree$tip.label)) %>%
     dplyr::select(phylum, colors) %>%
    unique() %>%
    arrange(phylum) %>%
    dplyr::select(colors) %>%
    pull()

# Basal ggtree
postestimates_tree <- postestimates_tree %>%
        force.ultrametric(.,method="extend") %>%
        ggtree(., size = 0.3)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
      scale_fill_manual(values=colors_alphabetic)+
      labs(fill="Phylum")

#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()

# Add posterior significant heatmap

postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
        scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
        labs(fill="Trend")

postestimates_tree +
        vexpand(.25, 1) # expand top 

8.3.1 Positively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Positive") %>%
    arrange(-value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum, class, family)%>%
    tt()
tinytable_hhg5ohm497z1vxt3adey
genome phylum class order family species value
bin_CaboVerde.50 Actinomycetota Actinomycetia Actinomycetales Actinomycetaceae NA 0.953
bin_Aruba.25 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium pseudocatenulatum 0.995
bin_Aruba.2 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium adolescentis 0.994
bin_Aruba.6 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium longum 0.985
bin_Denmark.14 Actinomycetota Actinomycetia Actinomycetales Bifidobacteriaceae Bifidobacterium gallinarum 0.950
bin_Aruba.47 Actinomycetota Actinomycetia Mycobacteriales Mycobacteriaceae NA 0.964
bin_Aruba.57 Actinomycetota Actinomycetia Mycobacteriales Mycobacteriaceae Corynebacterium pyruviciproducens 0.931
bin_Aruba.51 Actinomycetota Coriobacteriia Coriobacteriales Atopobiaceae Parolsenella sp900544655 0.969
bin_Denmark.44 Actinomycetota Coriobacteriia Coriobacteriales Atopobiaceae UBA7748 sp900314535 0.951
bin_Aruba.14 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 1.000
bin_Denmark.4 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella stercoris 1.000
bin_Aruba.21 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae NA 0.999
bin_Brazil.61 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella tanakaei 0.998
bin_Aruba.39 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp902362275 0.992
bin_Spain.84 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp900555555 0.992
bin_Brazil.151 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella intestinalis 0.991
bin_Denmark.127 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella sp900548365 0.991
bin_Denmark.17 Actinomycetota Coriobacteriia Coriobacteriales Coriobacteriaceae Collinsella bouchesdurhonensis 0.983
bin_Brazil.81 Bacillota Bacilli RFN20 CAG-826 NA 0.907
bin_Brazil.109 Bacillota_A Clostridia Oscillospirales Butyricicoccaceae AM07-15 sp003477405 0.936
bin_Aruba.28 Bacillota_A Clostridia Oscillospirales Oscillospiraceae NA 0.995
bin_Malaysia.9 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Dysosmobacter welbionis 0.986
bin_Malaysia.116 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Flavonifractor plautii 0.981
bin_Malaysia.34 Bacillota_A Clostridia Oscillospirales Oscillospiraceae Lawsonibacter sp000177015 0.954
bin_Aruba.54 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.959
bin_CaboVerde.35 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.930
bin_Aruba.36 Bacillota_A Clostridia Tissierellales Peptoniphilaceae Peptoniphilus_C sp902363535 0.913
bin_Aruba.52 Bacillota_A Clostridia Tissierellales Peptoniphilaceae NA 0.906
bin_Malaysia.26 Bacillota_A Clostridia Oscillospirales Ruminococcaceae NA 0.976
bin_Aruba.29 Bacillota_A Clostridia Oscillospirales Ruminococcaceae Gemmiger variabilis_C 0.941
bin_Malaysia.103 Bacillota_C Negativicutes Acidaminococcales Acidaminococcaceae Acidaminococcus intestini 0.942
bin_Malaysia.8 Bacillota_C Negativicutes Veillonellales Dialisteraceae Dialister sp900543165 0.981
bin_Denmark.60 Bacillota_C Negativicutes Veillonellales Dialisteraceae Allisonella histaminiformans 0.938
bin_Aruba.64 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera sp000417505 0.970
bin_Brazil.58 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera elsdenii 0.969
bin_CaboVerde.38 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera sp000417505 0.961
bin_Malaysia.81 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Megasphaera stantonii 0.943
bin_Malaysia.96 Bacillota_C Negativicutes Veillonellales Megasphaeraceae Caecibacter sp003467125 0.919
bin_Brazil.163 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella lascolaii 0.992
bin_CaboVerde.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella copri 0.981
bin_Aruba.10 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900544825 0.978
bin_Brazil.5 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotella sp900540415 0.929
bin_Spain.21 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola coprophilus 0.929
bin_Malaysia.18 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp000437675 0.924
bin_Malaysia.117 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Prevotellamassilia sp900541335 0.919
bin_Malaysia.77 Bacteroidota Bacteroidia Bacteroidales Bacteroidaceae Phocaeicola sp900542985 0.913
bin_Brazil.75 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae NA 0.964
bin_Malaysia.61 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae Helicobacter_C labetoulli 0.964
bin_Brazil.128 Campylobacterota Campylobacteria Campylobacterales Helicobacteraceae NA 0.963
bin_Brazil.70 Desulfobacterota Desulfovibrionia Desulfovibrionales Desulfovibrionaceae Mailhella sp003150275 0.944
bin_Brazil.146 Pseudomonadota Alphaproteobacteria RF32 CAG-239 CAG-495 sp001917125 0.969
bin_Brazil.9 Pseudomonadota Alphaproteobacteria RF32 CAG-239 CAG-495 sp000436375 0.919
bin_Aruba.15 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.963
bin_Brazil.51 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae Sutterella wadsworthensis_A 0.963
bin_Spain.43 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae CAG-521 sp000437635 0.959
bin_Brazil.64 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.948
bin_Brazil.92 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae CAG-521 sp000437635 0.948
bin_Brazil.15 Pseudomonadota Gammaproteobacteria Burkholderiales Burkholderiaceae NA 0.930
bin_CaboVerde.33 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum sp900543125 0.995
bin_Brazil.82 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum succiniciproducens 0.990
bin_CaboVerde.55 Pseudomonadota Gammaproteobacteria Enterobacterales Succinivibrionaceae Anaerobiospirillum_A thomasii 0.939

8.3.2 Negatively associated genomes

getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
    mutate(trend = case_when(
          value >= support ~ "Positive",
          value <= negsupport ~ "Negative",
          TRUE ~ "Neutral")) %>%
    filter(variable=="originTame", trend=="Negative") %>%
    arrange(value) %>%
    left_join(genome_metadata,by=join_by(genome==genome)) %>%
    dplyr::select(genome,phylum,class,order,family,species,value) %>%
    arrange(phylum,class,family)%>%
    tt()
tinytable_rk8eeroiunhr22x4voky
genome phylum class order family species value
bin_Denmark.27 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus_B hirae 0.000
bin_Brazil.12 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus faecalis 0.001
bin_Brazil.170 Bacillota Bacilli Lactobacillales Enterococcaceae NA 0.014
bin_CaboVerde.24 Bacillota Bacilli Lactobacillales Enterococcaceae Enterococcus_E cecorum 0.020
bin_CaboVerde.16 Bacillota Bacilli Lactobacillales Lactobacillaceae Limosilactobacillus reuteri 0.000
bin_Denmark.56 Bacillota Bacilli Lactobacillales Lactobacillaceae Levilactobacillus brevis 0.000
bin_Malaysia.72 Bacillota Bacilli Lactobacillales Lactobacillaceae Limosilactobacillus reuteri 0.000
bin_CaboVerde.60 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus agilis 0.001
bin_Malaysia.127 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus agilis 0.001
bin_Denmark.61 Bacillota Bacilli Lactobacillales Lactobacillaceae Latilactobacillus sakei 0.003
bin_Malaysia.35 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus animalis 0.003
bin_Malaysia.4 Bacillota Bacilli Lactobacillales Lactobacillaceae Ligilactobacillus ruminis 0.035
bin_CaboVerde.14 Bacillota Bacilli Lactobacillales Lactobacillaceae Lactobacillus acidophilus 0.067
bin_Aruba.18 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus lutetiensis 0.001
bin_Brazil.22 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus pasteurianus 0.001
bin_Denmark.117 Bacillota Bacilli Lactobacillales Streptococcaceae Lactococcus garvieae 0.001
bin_Denmark.113 Bacillota Bacilli Lactobacillales Streptococcaceae Streptococcus alactolyticus 0.002
bin_Brazil.69 Bacillota_A Clostridia Oscillospirales Acutalibacteraceae Clostridium_A leptum 0.097
bin_Denmark.93 Bacillota_A Clostridia Lachnospirales Anaerotignaceae Anaerotignum sp001304995 0.074
bin_CaboVerde.3 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium_P perfringens 0.000
bin_Denmark.42 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium_P perfringens 0.000
bin_Brazil.136 Bacillota_A Clostridia Clostridiales Clostridiaceae Clostridium thermobutyricum 0.041
bin_Brazil.3 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900555395 0.026
bin_Brazil.89 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900550235 0.034
bin_Denmark.63 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerostipes hadrus 0.040
bin_Denmark.19 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Faecalimonas sp900556835 0.052
bin_Denmark.118 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Dorea_A longicatena 0.067
bin_Malaysia.108 Bacillota_A Clostridia Lachnospirales Lachnospiraceae NA 0.081
bin_Brazil.113 Bacillota_A Clostridia Lachnospirales Lachnospiraceae Anaerobutyricum sp900754855 0.095
bin_Denmark.50 Bacillota_A Clostridia Peptostreptococcales Peptostreptococcaceae Peptacetobacter sp900539645 0.053
bin_Brazil.107 Bacillota_A Clostridia Monoglobales_A UBA1381 CAG-41 sp900066215 0.016
bin_Brazil.159 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900541465 0.040
bin_Brazil.140 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_A sp900543175 0.046
bin_Malaysia.6 Fusobacteriota Fusobacteriia Fusobacteriales Fusobacteriaceae Fusobacterium_B sp900545035 0.088

8.3.3 Estimate - support plot

estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
    dplyr::select(genome,mean)

support <- getPostEstimate(hM=m, parName="Beta")$support %>%
    as.data.frame() %>%
    mutate(variable=m$covNames) %>%
    filter(variable=="originTame") %>%
    pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
    dplyr::select(genome,support)

phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
  right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
  filter(genome %in% estimate$genome) %>%
  arrange(match(genome, estimate$genome)) %>%
  dplyr::select(phylum, colors) %>%
  unique() %>%
  arrange(phylum) %>%
  dplyr::select(colors) %>%
  pull()
Rows: 202 Columns: 2
── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): phylum, colors

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
inner_join(estimate,support,by=join_by(genome==genome)) %>%
    mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
    mutate(support=ifelse(mean<0,1-support,support)) %>%
    left_join(genome_metadata, by = join_by(genome == genome)) %>%
    mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
    ggplot(aes(x=mean,y=support,color=phylum))+
      geom_point(alpha=0.7, shape=16, size=3)+
      scale_color_manual(values = phylum_colors) +
      geom_vline(xintercept = 0) +
      xlim(c(-0.4,0.4)) +
      labs(x = "Beta regression coefficient", y = "Posterior probability") +
      theme_minimal()+
       theme(legend.position = "none")

8.4 Correlations

#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)

# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
        keep.tip(., tip=m$spNames) 


#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
          + (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean

matrix <- toPlot %>% 
      as.data.frame() %>%
      rownames_to_column(var="genome1") %>%
      pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
      mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
      mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
      ggplot(aes(x = genome1, y = genome2, fill = cor)) +
            geom_tile() + 
            scale_fill_gradient2(low = "#be3e2b",
                       mid = "#f4f4f4",
                       high = "#b2b530")+
            theme_void()

htree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
vtree <- genome_tree_subset %>%
  force.ultrametric(.,method="extend") %>%
  ggtree(.)
***************************************************************
*                          Note:                              *
*    force.ultrametric does not include a formal method to    *
*    ultrametricize a tree & should only be used to coerce    *
*   a phylogeny that fails is.ultrametric due to rounding --  *
*    not as a substitute for formal rate-smoothing methods.   *
***************************************************************
#create composite figure
grid.arrange(grobs = list(matrix,vtree),
             layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1),
                                   c(2,1,1,1,1,1,1,1,1,1,1,1)))

8.5 Predict responses

# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")

gradient = c("domestic","feral")
gradientlength = length(gradient)

#Treatment-specific gradient predictions
pred <- constructGradient(m, 
                      focalVariable = "origin", 
                      non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
            predict(m, Gradient = ., expected = TRUE) %>%
            do.call(rbind,.) %>%
            as.data.frame() %>%
            mutate(origin=rep(gradient,1000)) %>%
            pivot_longer(!origin,names_to = "genome", values_to = "value")
# weights:  9 (4 variable)
initial  value 101.072331 
final  value 91.392443 
converged

8.5.0.1 Element level

elements_table <- genome_gifts %>%
    to.elements(., GIFT_db) %>%
    as.data.frame()

community_elements <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(elements_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })

calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

element_predictions <- map_dfc(community_elements, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_elements[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

8.5.1 Positive associated functions at element level

# Positively associated

unique_funct_db<- GIFT_db[c(2,4,5,6)] %>% 
  distinct(Code_element, .keep_all = TRUE)


element_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  filter(positive_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_ytjmkzr7zycoaze18ev3
trait mean p1 p9 positive_support negative_support Domain Function Element
B0103 0.008721055 5.858948e-05 0.016907646 0.902 0.098 Biosynthesis Nucleic acid biosynthesis UDP/UTP
D0504 0.004386196 1.313410e-04 0.009500262 0.911 0.089 Degradation Amino acid degradation Methionine
D0906 0.003670175 1.379577e-04 0.008121394 0.926 0.074 Degradation Antibiotic degradation Fosfomycin
D0205 0.012779324 2.341355e-03 0.023302092 0.954 0.046 Degradation Polysaccharide degradation Pectin
D0208 0.010331977 1.709235e-03 0.018996311 0.935 0.065 Degradation Polysaccharide degradation Mixed-Linkage glucans
element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_mssyxcyyf5ercqt95n2s
trait mean p1 p9 positive_support negative_support Domain Function Element
B0219 -0.004243274 -0.009449620 -2.360171e-04 0.041 0.959 Biosynthesis Amino acid biosynthesis GABA
B0214 -0.020976378 -0.036675942 -4.271212e-03 0.069 0.931 Biosynthesis Amino acid biosynthesis Glutamate
B0302 -0.005039082 -0.010242469 -5.152291e-04 0.032 0.968 Biosynthesis Amino acid derivative biosynthesis Betaine
B0310 -0.012921300 -0.022832185 -3.166348e-03 0.037 0.963 Biosynthesis Amino acid derivative biosynthesis Tryptamine
B0303 -0.011717187 -0.021302206 -2.302478e-03 0.059 0.941 Biosynthesis Amino acid derivative biosynthesis Ectoine
B0309 -0.007622324 -0.014859566 -3.193213e-04 0.092 0.908 Biosynthesis Amino acid derivative biosynthesis Putrescine
B0804 -0.016338997 -0.029159443 -2.639667e-03 0.057 0.943 Biosynthesis Aromatic compound biosynthesis Dipicolinate
B0603 -0.016187723 -0.031768524 -2.015090e-03 0.076 0.924 Biosynthesis Organic anion biosynthesis Citrate
B0601 -0.009214752 -0.018849619 -8.066553e-04 0.080 0.920 Biosynthesis Organic anion biosynthesis Succinate
B0401 -0.011952355 -0.022575997 -1.074075e-03 0.071 0.929 Biosynthesis SCFA biosynthesis Acetate
B0709 -0.002065553 -0.003638992 -5.805158e-04 0.032 0.968 Biosynthesis Vitamin biosynthesis Tocopherol/tocotorienol (E)
B0708 -0.024292457 -0.048067418 -2.361969e-04 0.098 0.902 Biosynthesis Vitamin biosynthesis Cobalamin (B12)
D0517 -0.004760845 -0.008589668 -1.196185e-03 0.031 0.969 Degradation Amino acid degradation Ornithine
D0508 -0.003427595 -0.007405302 -2.064998e-04 0.074 0.926 Degradation Amino acid degradation Lysine
D0903 -0.004207352 -0.009395988 -2.301420e-04 0.041 0.959 Degradation Antibiotic degradation Cephalosporin
D0908 -0.015747946 -0.028054674 -3.756246e-03 0.068 0.932 Degradation Antibiotic degradation Macrolide
D0601 -0.009353606 -0.017995591 -2.103861e-03 0.034 0.966 Degradation Nitrogen compound degradation Nitrate
D0611 -0.004207352 -0.009395988 -2.301420e-04 0.041 0.959 Degradation Nitrogen compound degradation Phenylethylamine
D0603 -0.001979236 -0.003939691 -3.288703e-04 0.043 0.957 Degradation Nitrogen compound degradation Urate
D0610 -0.002989440 -0.004961091 -8.944614e-04 0.052 0.948 Degradation Nitrogen compound degradation Methylamine
D0606 -0.005920775 -0.011186099 -1.007489e-03 0.061 0.939 Degradation Nitrogen compound degradation Allantoin
D0612 -0.001692843 -0.003140507 -1.717185e-05 0.097 0.903 Degradation Nitrogen compound degradation Hypotaurine
D0801 -0.001729950 -0.002274583 -9.185503e-05 0.012 0.988 Degradation Xenobiotic degradation Toluene
D0802 -0.001729950 -0.002274583 -9.185503e-05 0.012 0.988 Degradation Xenobiotic degradation Xylene
D0807 -0.004248376 -0.008547901 -6.633122e-04 0.054 0.946 Degradation Xenobiotic degradation Catechol
D0817 -0.005084664 -0.010358431 -5.854535e-04 0.059 0.941 Degradation Xenobiotic degradation Trans-cinnamate
D0816 -0.005861200 -0.012171134 -8.084910e-05 0.093 0.907 Degradation Xenobiotic degradation Phenylacetate

8.5.2 Negatively associated functions at element level

element_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  filter(negative_support>=0.9) %>%
  left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
  arrange(Domain,Function)%>%
  tt()
tinytable_ro9ms6tz0qt2gx9au5su
trait mean p1 p9 positive_support negative_support Domain Function Element
B0219 -0.004243274 -0.009449620 -2.360171e-04 0.041 0.959 Biosynthesis Amino acid biosynthesis GABA
B0214 -0.020976378 -0.036675942 -4.271212e-03 0.069 0.931 Biosynthesis Amino acid biosynthesis Glutamate
B0302 -0.005039082 -0.010242469 -5.152291e-04 0.032 0.968 Biosynthesis Amino acid derivative biosynthesis Betaine
B0310 -0.012921300 -0.022832185 -3.166348e-03 0.037 0.963 Biosynthesis Amino acid derivative biosynthesis Tryptamine
B0303 -0.011717187 -0.021302206 -2.302478e-03 0.059 0.941 Biosynthesis Amino acid derivative biosynthesis Ectoine
B0309 -0.007622324 -0.014859566 -3.193213e-04 0.092 0.908 Biosynthesis Amino acid derivative biosynthesis Putrescine
B0804 -0.016338997 -0.029159443 -2.639667e-03 0.057 0.943 Biosynthesis Aromatic compound biosynthesis Dipicolinate
B0603 -0.016187723 -0.031768524 -2.015090e-03 0.076 0.924 Biosynthesis Organic anion biosynthesis Citrate
B0601 -0.009214752 -0.018849619 -8.066553e-04 0.080 0.920 Biosynthesis Organic anion biosynthesis Succinate
B0401 -0.011952355 -0.022575997 -1.074075e-03 0.071 0.929 Biosynthesis SCFA biosynthesis Acetate
B0709 -0.002065553 -0.003638992 -5.805158e-04 0.032 0.968 Biosynthesis Vitamin biosynthesis Tocopherol/tocotorienol (E)
B0708 -0.024292457 -0.048067418 -2.361969e-04 0.098 0.902 Biosynthesis Vitamin biosynthesis Cobalamin (B12)
D0517 -0.004760845 -0.008589668 -1.196185e-03 0.031 0.969 Degradation Amino acid degradation Ornithine
D0508 -0.003427595 -0.007405302 -2.064998e-04 0.074 0.926 Degradation Amino acid degradation Lysine
D0903 -0.004207352 -0.009395988 -2.301420e-04 0.041 0.959 Degradation Antibiotic degradation Cephalosporin
D0908 -0.015747946 -0.028054674 -3.756246e-03 0.068 0.932 Degradation Antibiotic degradation Macrolide
D0601 -0.009353606 -0.017995591 -2.103861e-03 0.034 0.966 Degradation Nitrogen compound degradation Nitrate
D0611 -0.004207352 -0.009395988 -2.301420e-04 0.041 0.959 Degradation Nitrogen compound degradation Phenylethylamine
D0603 -0.001979236 -0.003939691 -3.288703e-04 0.043 0.957 Degradation Nitrogen compound degradation Urate
D0610 -0.002989440 -0.004961091 -8.944614e-04 0.052 0.948 Degradation Nitrogen compound degradation Methylamine
D0606 -0.005920775 -0.011186099 -1.007489e-03 0.061 0.939 Degradation Nitrogen compound degradation Allantoin
D0612 -0.001692843 -0.003140507 -1.717185e-05 0.097 0.903 Degradation Nitrogen compound degradation Hypotaurine
D0801 -0.001729950 -0.002274583 -9.185503e-05 0.012 0.988 Degradation Xenobiotic degradation Toluene
D0802 -0.001729950 -0.002274583 -9.185503e-05 0.012 0.988 Degradation Xenobiotic degradation Xylene
D0807 -0.004248376 -0.008547901 -6.633122e-04 0.054 0.946 Degradation Xenobiotic degradation Catechol
D0817 -0.005084664 -0.010358431 -5.854535e-04 0.059 0.941 Degradation Xenobiotic degradation Trans-cinnamate
D0816 -0.005861200 -0.012171134 -8.084910e-05 0.093 0.907 Degradation Xenobiotic degradation Phenylacetate
positive <- element_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- element_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

table <- bind_rows(positive,negative) %>%
  left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) 

table %>%
  mutate(Element=factor(Element,levels=c(table$Element))) %>%
  ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.04,0.04)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")

8.5.2.1 Function level

functions_table <- elements_table %>%
    to.functions(., GIFT_db) %>%
    as.data.frame()

community_functions <- pred %>%
  group_by(origin, genome) %>%
  mutate(row_id = row_number()) %>%
  pivot_wider(names_from = genome, values_from = value) %>%
  ungroup() %>%
  group_split(row_id) %>%
  as.list() %>%
  lapply(., FUN = function(x){x %>%
    dplyr::select(-row_id) %>%
    column_to_rownames(var = "origin") %>%
    as.data.frame() %>%
    exp() %>%
    t() %>%
    tss() %>%
    to.community(functions_table,.,GIFT_db) %>% 
    as.data.frame() %>%
    rownames_to_column(var="origin")
   })
#max-min option
calculate_slope <- function(x) {
  lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
  coef(lm_fit)[2]
}

function_predictions <- map_dfc(community_functions, function(mat) {
      mat %>%
        column_to_rownames(var = "origin") %>%
        t() %>%
        as.data.frame() %>%
        rowwise() %>%
        mutate(slope = calculate_slope(c_across(everything()))) %>%
        dplyr::select(slope) }) %>%
      t() %>%
      as.data.frame() %>%
      set_names(colnames(community_functions[[1]])[-1]) %>%
      rownames_to_column(var="iteration") %>%
      pivot_longer(!iteration, names_to="trait",values_to="value") %>%
      group_by(trait) %>%
      summarise(mean=mean(value),
        p1 = quantile(value, probs = 0.1),
        p9 = quantile(value, probs = 0.9),
        positive_support = sum(value > 0)/1000,
        negative_support = sum(value < 0)/1000) %>%
      arrange(-positive_support)

# Positively associated
function_predictions %>%
  filter(mean >0) %>%
  arrange(-positive_support) %>%
  tt()
tinytable_d4yj5jt63t55yrd0cyr1
trait mean p1 p9 positive_support negative_support
D02 8.580932e-03 -0.0029113778 0.0211120063 0.826 0.174
B08 7.733393e-03 -0.0033231015 0.0179292961 0.792 0.208
B01 1.376092e-03 -0.0058174809 0.0084040031 0.647 0.353
S01 8.140471e-04 -0.0125251212 0.0146119095 0.572 0.428
B09 5.381937e-05 -0.0005426244 0.0005120737 0.371 0.629
# Negatively associated
function_predictions %>%
  filter(mean <0) %>%
  arrange(-negative_support) %>%
  tt()
tinytable_g38s5y613luak78vxnn0
trait mean p1 p9 positive_support negative_support
D08 -1.174084e-03 -0.0021635400 -0.0002216101 0.047 0.953
B03 -1.082132e-02 -0.0184920050 -0.0036627348 0.054 0.946
D06 -3.194052e-03 -0.0070018415 0.0001290343 0.116 0.884
B04 -8.373734e-03 -0.0174614803 0.0007363593 0.121 0.879
D07 -1.200153e-02 -0.0300385854 0.0042494388 0.185 0.815
B06 -6.669703e-03 -0.0166283189 0.0027105789 0.199 0.801
D05 -1.717352e-03 -0.0075185682 0.0038509242 0.216 0.784
S03 -9.563521e-03 -0.0333823692 0.0159823031 0.226 0.774
D03 -4.017320e-03 -0.0127303011 0.0040172389 0.236 0.764
B02 -3.462674e-03 -0.0124506768 0.0051573014 0.298 0.702
D09 -2.167099e-03 -0.0085278347 0.0046009115 0.299 0.701
S02 -4.696911e-03 -0.0151796566 0.0030908687 0.323 0.677
B07 -3.645464e-03 -0.0169431855 0.0096020470 0.361 0.639
D01 -3.640355e-04 -0.0053797539 0.0041408263 0.422 0.578
B10 -3.069883e-06 -0.0002900964 0.0002734412 0.473 0.527
positive <- function_predictions %>%
  filter(mean >0) %>%
  arrange(mean) %>%
  filter(positive_support>=0.9) %>%
  dplyr::select(-negative_support) %>%
  rename(support=positive_support)

negative <- function_predictions %>%
  filter(mean <0) %>%
  arrange(mean) %>%
  filter(negative_support>=0.9) %>%
  dplyr::select(-positive_support) %>%
  rename(support=negative_support)

bind_rows(positive,negative) %>%
  left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
  mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
  ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
      geom_point() +
      geom_errorbar() +
      xlim(c(-0.02,0.02)) +
      geom_vline(xintercept=0) +
      scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
      theme_minimal() +
      labs(x="Regression coefficient",y="Functional trait")